

package scalaSci

import scalaSci.math.LinearAlgebra.LinearAlgebra
import scalaSci.math.array.DoubleArray
import java.util._
import Jama._


// this class provides Matrix operations, uses one indexed values for direct interfacing with NumAL
class Matrix( n: Int, m: Int) {    // the default constructor allocates a double array of size n x m
   private var  N: Int=n+1
   private var  M: Int=m+1
   private  var  v = new Array[Array[Double]](N, M)  // i.e. v(1:n, 1:m)

  def  setv(values: Array[Array[Double]], n: Int, m: Int) = { v = values; N = n+1; M = m+1 }
   
  def size = {  (N-1, M-1) }  // the size of the Matrix

  def length =  { N-1 }

  // getters for size
  def Nrows = { N-1 }
  def Ncols = { M-1 }

  def getv =    v

    // construct a Matrix from a zero-indexed double [][] array
def this(vals: Array[Array[Double]]) = {
     this(vals.length, vals(0).length)  // allocate memory with the default constructor
    var i=0; var j=0
    while (i<vals.length) {
        j=0
        while (j<vals(0).length) {
       v(i+1)(j+1) = vals(i)(j)
       j+=1
        }
        i+=1
     }
 }



def this(N: Int, M: Int, df: Double) = {   // initializes to a default value
     this(N, M)  // allocate memory with the default constructor
    var i=0; var j=0
    while (i<v.length-1) {
        j=0
        while (j<v(1).length-1) {
       v(i+1)(j+1) = df
       j+=1
        }
        i+=1
     }
 }

    // construct a Matrix from an one-indexed double [][] array
def this(vals: Array[Array[Double]], flag: Boolean) = {
    this(1,1)
    v = vals
    N = vals.length
    M = vals(0).length
 }

  // construct a Matrix from a one-indexed double [] array
def this( vals: Array[Double]) = {
    this(1, vals.length)   // keep the array as a row of the Matrix
    var i=1
    while (i<vals.length)  {
       v(1)(i) = vals(i)
       i+=1
   }
  }

  // construct a Matrix from a double [] array
def this( vals: Array[Double], oneIndexed: Boolean) = {
    this(1, vals.length-1)   // keep the array as a row of the Matrix 
    var i=1
    while (i<vals.length)  {
       v(1)(i) = vals(i)
       i+=1
   }
  }

// apply the function f to all the elements of the Matrix and return the results with a new Matrix
 def  map( f: (Double => Double)): Matrix = {
   var mres = new Matrix(N, M)
   
    for (r<-0 until N)
      for (c<-0 until M) 
       mres(r, c) = f(this(r, c) )
   
   mres
 }

  // clone a Matrix from another Matrix
def this(a: Matrix) = {
     this(a.N-1, a.M-1)  // allocate memory with the default constructor
    var i=1; var j=1
    while (i<a.N) {
        j=1
        while (j<a.M) {
       v(i)(j) = a.v(i)(j)
       j+=1
        }
        i+=1
     }
 }
     
def  print: String = {
  var digitFormat = scalaExec.Interpreter.GlobalValues.fmtMatrix
    var nelems = N
    var melems = M

    var s:String =" \n"
     var i=1; var j=1;
    while (i<nelems) {
        j=1
        while (j<melems) {
       s = s+digitFormat.format(v(i)(j))+"  "
       j+=1
        }
        s += "\n";
     i+=1
    }
  s
 }

def p = print   // short "print" method


override  def  toString(): String = {
  var digitFormat = scalaExec.Interpreter.GlobalValues.fmtMatrix
    var nelems = N
    var melems = M
    var truncated = false

    var truncationIndicator = ""
     if  (Matrix.mxElemsToDisplay < N) {
        nelems = Matrix.mxElemsToDisplay
        truncationIndicator = " ... "
        truncated = true
      }
     if  (Matrix.mxElemsToDisplay < M) {
        melems = Matrix.mxElemsToDisplay
        truncationIndicator = " ... "
      }
     
    var s:String =" \n"
     var i=1; var j=1;
    while (i<nelems) {
        j=1
        while (j<melems) {
       s = s+digitFormat.format(v(i)(j))+"  "
       j+=1
        }
     s += (truncationIndicator+"\n")
     i+=1
    }
 if (truncated)     // handle last line
    s+= "........................................................................................................................................."
     return s
 }
  

// indexes the corresponding Matrix element
  def apply(n: Int, m: Int): Double = {
          v(n)(m)
 }

 
  // indexes the corresponding Matrix element, resizing automatically 
  def apply(n: Int, m: Int, resize: Boolean): Double = {
      if (n>=N || m >= M)  {  // update size
          var newN = (Matrix.matResizeFactor*n).asInstanceOf[Int]
          var newM = (Matrix.matResizeFactor*m).asInstanceOf[Int]
              // create the new larger matrix
          var v2 = new Array[Array[Double]](newN)
          var c=0
          while (c<newN) {
              v2(c) = new Array[Double](newM)
              c += 1
          }
            // copy the old one
            var i=0;  var j=0
            while (i<N)  {
                j=0;
                while (j<M)   {
                    v2(i)(j) = v(i)(j)
                    j += 1
                  }
                i += 1
            }
            // update the matrix representation to  the larger one
          v = v2
          N = newN
          M = newM
      }   // update size

        v(n)(m)

 }

 
// extracts a submatrix specifying rows only, take all columns, e.g. m(2, 3, ':') corresponds to Matlab's m(2:3, :)'
// m(low:high,:) is implemented with m(low, high, dummyChar). if low>high then rows are returned in reverse
 def apply(rowL: Int, rowH: Int, allColsChar: Char): Matrix = {
   var rowStart = rowL; var rowEnd=rowH;
   var colStart = 1;     var colEnd =  M-1;   // all columns
   var colNum = M-1
   var colInc = 1

if (rowStart <= rowEnd) {   // positive increment
    var rowInc = 1   
    if (rowEnd == -1) { rowEnd = N-1 }  // if -1 is specified take all the rows
    var rowNum = rowEnd-rowStart+1 
    var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart
    var rowIdx =1; var colIdx = 1  // indexes at the new Matrix
    while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol <= colEnd )   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
       } // crow <= rowEnd
subMatr   // return the submatrix
        
} // rowStart <= rowEnd
else { // rowStart > rowEnd
    var rowInc = -1   
    var rowNum = rowStart-rowEnd+1
    var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row at the source matrix
    var ccol = colStart
    var rowIdx =1; var colIdx = 1  // indexes at the new Matrix
    while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
       }

subMatr   // return the submatrix
        
} // rowStart > rowEnd
  
}

// extracts a submatrix specifying rows only, take all columns, e.g. m(2, 4, 12, ':') corresponds to Matlab's m(2:4:12, :)'
def apply(rowL: Int, rowInc: Int, rowH: Int, allColsChar: Char): Matrix = {
    var rowStart = rowL;     var rowEnd =  rowH;
    var colStart = 1;  var colEnd = M-1;   // all columns
    var colNum = M-1
    var colInc = 1

  if (rowInc > 0) { // positive increment
    if (rowEnd == -1) { rowEnd = N-1 }  // if -1 is specified take all the rows
    var rowNum = Math.floor( (rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
    var colStart = 1;     var colEnd =  M-1   // all columns
    var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix
    while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
       }
     subMatr   // return the submatrix
     }  // positive increment
  else  {  //  negative increment
     var rowNum = Math.floor( (rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
     var subMatr = new Matrix(rowNum, colNum)  // create a Matrix to keep the extracted range
        // fill the created matrix with values
     var  crow = rowStart   // indexes current row at the source matrix
     var  ccol = colStart
     var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix
     while (crow >= rowEnd)  {
         ccol = colStart;  colIdx = 1
         while (ccol <= colEnd)  {
             subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
             colIdx += 1
             ccol += colInc
         }
         rowIdx += 1
         crow += rowInc
      }
       subMatr  // return the submatrix
     }  // negative increment
}


// extracts a submatrix, e.g. m( ':', 2, 3 ) corresponds to Matlab's m(:, 2:3)'
  def apply(allRowsChar: Char, colLow: Int, colHigh: Int): Matrix = {
    var rowStart = 1;     var rowEnd =  N-1   // all rows
    var colStart = colLow;  var colEnd = colHigh
    var rowInc=1
    var rowNum = N-1    // take all the rows

    if  (colStart <= colEnd)   {    // positive increment
        var colInc = 1
        if (colEnd == -1)  { colEnd = M-1 } // if -1 is specified take all the columns
        var colNum = colEnd-colStart+1;
    var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEbd
 subMatr
} // positive increment
  else {  // negative increment
    var colInc = -1
    var colNum = colStart-colEnd+1;
    var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
   
  }



// extracts a submatrix, e.g. m( ':', 2, 3, 12 ) corresponds to Matlab's m(:, 2:3:12)'
  def apply(allRowsChar: Char, colLow: Int, colInc: Int, colHigh: Int): Matrix = {
   var rowStart = 1;     var rowEnd =  N-1   // all rows
    var colStart = colLow;  var colEnd = colHigh
    var rowInc=1
    var rowNum = N-1    // take all the rows

    if  (colStart <= colEnd)   {    // positive increment
        if (colEnd == -1)  { colEnd = M-1 } // if -1 is specified take all the columns
        var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
        var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix
 
           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr
} // positive increment
  else {  // negative increment
    var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
    var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
   }


// extracts a submatrix, e.g. m( 2, 3, 12, 4, 2,  8 ) corresponds to Matlab's m(2:3:12, 4:2:8)'
  def apply(rowLow: Int, rowInc: Int, rowHigh: Int, colLow: Int, colInc: Int, colHigh: Int): Matrix = {
    var rowStart = rowLow;     var rowEnd =  rowHigh
    var colStart = colLow;  var colEnd = colHigh
    
        var rowNum = Math.floor((rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
        var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
        var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      
    if  (rowStart <= rowEnd && colStart <= colEnd)   {    // positive increment at rows and columns
        var crow = rowStart  // indexes current row
        var ccol = colStart  // indexes current column
        var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix
            while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr
} // positive increment
  else if  (rowStart >= rowEnd && colStart <= colEnd)   {    
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix

           while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
else if  (rowStart <= rowEnd && colStart >= colEnd)   {    
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
else {
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix

           while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow > rowEnd
 subMatr   // return the submatrix
   }

   }



// extracts a specific row, take all columns, e.g. m(2) corresponds to Matlab's m(2, :)'
  def apply(row: Int): Matrix = {
    var colStart = 1;     var colEnd =  M-1;   // all columns
    var rowNum = 1;  var colNum = colEnd-colStart+1;
    var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var ccol = colStart
    while  (ccol <= colEnd)   {
          subMatr.v(1)(ccol) = v(row)(ccol)
          ccol += 1
         }
       
     subMatr
}


// extracts a specific column, take all rows, e.g. m(':', 2) corresponds to Matlab's m(:,2:)'
  def apply(colonChar:Char, col: Int): Matrix = {
    var rowStart = 1;     var rowEnd =  N-1;   // all rows
    var colNum = 1;  var rowNum = rowEnd-rowStart+1;
    var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart
    while  (crow <= rowEnd)   {
          subMatr.v(crow)(1) = v(crow)(col)
          crow += 1
         }
       
     subMatr
}

// extracts a specific row, take all columns, e.g. m(2, ':') corresponds to Matlab's m(2, :)'
  def apply(row: Int, allColsChar: Char): Matrix = {
      if (allColsChar==':')
          apply(row)
      this
      }


   
// extracts a submatrix, e.g. m( 2,  12, 4,   8 ) corresponds to Matlab's m(2:12, 4:8)'
  def apply(rowLow: Int,  rowHigh: Int, colLow: Int, colHigh: Int): Matrix = {
    var rowStart = rowLow;     var rowEnd =  rowHigh
    var colStart = colLow;  var colEnd = colHigh
    var rowInc = if (rowHigh > rowLow) 1 else -1
    var colInc = if (colHigh > colLow) 1 else -1

        var rowNum = Math.floor((rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
        var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
        var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      
    if  (rowStart <= rowEnd && colStart <= colEnd)   {    // positive increment at rows and columns
        var crow = rowStart  // indexes current row
        var ccol = colStart  // indexes current column
        var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix
            while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr
} // positive increment
  else if  (rowStart >= rowEnd && colStart <= colEnd)   {    
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix

           while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol <= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
else if  (rowStart <= rowEnd && colStart >= colEnd)   {    
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix

           while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow <= rowEnd
 subMatr   // return the submatrix
   }
else {
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var ccol = colStart  // indexes current column
    var rowIdx = 1; var colIdx = 1  // indexes at the new Matrix

           while  ( crow >= rowEnd )   {
          ccol = colStart;  colIdx = 1
          while  (ccol >= colEnd)   {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
     } // crow > rowEnd
 subMatr   // return the submatrix
   }

   }




              // extracts a submatrix, e.g. m(3:2:7, :)
  def apply(rowLow: Int, rowInc: Int, rowHigh: Int): Matrix = {
    var rowStart = rowLow;     var rowEnd =  rowHigh;    if (rowEnd < rowStart) { rowStart = rowHigh; rowEnd = rowLow; }
    var colStart = 1;     var colEnd =  M-1;
    var colInc = 1
    var rowNum = Math.floor( (rowEnd-rowStart) / rowInc).asInstanceOf[Int]+1
    var colNum = Math.floor( (colEnd-colStart) / colInc).asInstanceOf[Int]+1
    var subMatr = new Matrix(rowNum, colNum)   // create a Matrix to keep the extracted range
      // fill the created matrix with values
    var crow = rowStart  // indexes current row
    var inc = 1
    var ccol = colStart
    var rowIdx =1; var colIdx = 1;  // indexes at the new Matrix
    while  ( crow <= rowEnd )   {
          ccol = colStart;  colIdx = 0;
          while  (ccol <= colEnd)    {
                subMatr.v(rowIdx)(colIdx) = v(crow)(ccol)
                colIdx += 1
                ccol += colInc
               }
            rowIdx += 1
            crow += rowInc
       }
     subMatr
}


// returns the corresponding row of the Matrix class as an Array[Double]
def getRow(row: Int): Array[Double] = {
    var colStart = 0;     var colEnd =  M-1;   // all columns
    var colNum = colEnd-colStart+1;
    var rowArray = new Array[Double](colNum)
    for (ccol<-0 to colEnd)
       rowArray(ccol) = v(row)(ccol)
    rowArray
  }


  // updating a single element of the Matrix, without resizing
def  update(n: Int, m: Int, value: Double ): Unit = {
         v(n)(m) = value;
   }

  
   // updating a single element of the Matrix, automatic resizing
def  update(n: Int, m: Int, value: Double, resize: Boolean ): Unit = {
   if (n>=N || m >= M)  {  // update size
          var newN = (Matrix.matResizeFactor*n).asInstanceOf[Int]
          var newM = (Matrix.matResizeFactor*m).asInstanceOf[Int]
              // create the new larger matrix
          var v2 = new Array[Array[Double]](newN)
          var c=1
          while (c<newN) {
              v2(c) = new Array[Double](newM)
              c += 1
          }
            // copy the old one
            var i=1;  var j=1
            while (i<N)  {
                j=1
                while (j<M)   {
                    v2(i)(j) = v(i)(j)
                    j += 1
                  }
                i += 1
            }
            // update the matrix representation to  the larger one
          v = v2
          N = newN
          M = newM
      }   // update size

         v(n)(m) = value;
   }

    // updating a single element of the Matrix. no resizing
def  update(n: Int, m: Int, value: Double, resizeDummy: Double): Unit = {
        v(n)(m) = value;
   }

// update a Matrix subrange by assigning a Matrix, e.g. var m=rand(20, 30);  m( 2, 1, 3, 3, 1, 4) = ones(2,2)
 def update(rlowp:Int, rincp: Int, rhighp: Int, clowp:Int, cincp: Int,  chighp: Int, mr: Matrix ): Unit = {
     var rlow=rlowp; var rhigh=rhighp; var rinc = rincp;
     var clow=clowp; var chigh=chighp; var cinc = cincp;
     if (rhigh < rlow)  {
            if (rinc > 0)  {
                println("negative row subrange increment is required")
                this
            }
            var tmp=rhighp; var rhigh=rlowp; rlow = tmp;
            rinc  = -rinc
        }
    if (chigh < clow)  {
            if (cinc > 0)  {
                println("negative column subrange increment is required")
                this
            }
            var tmp=chighp; var chigh=clowp; clow = tmp;
            cinc  = -cinc
        }
      var mrM = mr.M-1   // length of right-hand side matrix
      var mrN = mr.N-1
     var rangeLenRow = ((rhigh-rlow)/rinc).asInstanceOf[Int]+1    // row length of target range
     var rangeLenCol = ((chigh-clow)/cinc).asInstanceOf[Int]+1    // col length of target range

     if (mrN != rangeLenRow)  {  // improper sizes of the involved row subranges for assignment
         println("target vector row subrange and right-hand size vector lengths do not agree, i.e. mrN =  "+mrN+", rangeLenRow = "+rangeLenRow)
         this   // return the vector unaltered
     }

     if (mrM != rangeLenCol)  {  // improper sizes of the involved column subranges for assignment
         println("target vector column subrange and right-hand size vector lengths do not agree, i.e. mrM =  "+mrM+", rangeLenRow = "+rangeLenRow)
         this   // return the vector unaltered
     }

      if (rhigh >= N || chigh >= M)  {   // dynamically increase the size of the matrix when subassigning out of its range
          var newN = (Matrix.matResizeFactor*rhigh).asInstanceOf[Int]
          var newM = (Matrix.matResizeFactor*chigh).asInstanceOf[Int]
              // create the new larger matrix
          var v2 = new Array[Array[Double]](newN)
          var c=0
          while (c<newN) {
              v2(c) = new Array[Double](newM)
              c += 1
          }
            // copy the old one
            var i=0;  var j=0
            while (i<N)  {
                j=0;
                while (j<M)   {
                    v2(i)(j) = v(i)(j)
                    j += 1
                  }
                i += 1
            }
            // update the matrix representation to  the larger one
          v = v2
          N = newN
          M = newM
      }   // dynamically increase the size of the matrix


     // copy the values of the mr
        var rrow=1; var rcol=1; var lrowidx=rlow; var lcolidx = clow
        while (rrow < mr.N) {   // for all rows of the right-hand side matrix
            rcol=1
            lcolidx = clow   // starting column within the "subassigned" matrix
            while (rcol < mr.M)  {   // for all cols of the right-hand side matrix
                v(lrowidx)(lcolidx) = mr(rrow, rcol)
                lcolidx += cinc
                rcol += 1
            }
            lrowidx += rinc
            rrow += 1
          }
            
        }



// update a Matrix subrange by assigning a Matrix, e.g. var mm = rand(20, 30);  mm(2,3, 3,4) = ones(2,2);
 def update(rlowp:Int, rhighp: Int, clowp:Int, chighp: Int, mr: Matrix ): Unit = {
     var rlow=rlowp; var rhigh=rhighp; var rinc = 1;
     var clow=clowp; var chigh=chighp; var cinc = 1;
     if (rhigh < rlow)  {
            if (rinc > 0)  {
                println("negative row subrange increment is required")
                this
            }
            var tmp=rhighp; var rhigh=rlowp; rlow = tmp;
            rinc  = -rinc
        }
    if (chigh < clow)  {
            if (cinc > 0)  {
                println("negative column subrange increment is required")
                this
            }
            var tmp=chighp; var chigh=clowp; clow = tmp;
            cinc  = -cinc
        }
      var mrM = mr.M   // length of right-hand side matrix
      var mrN = mr.N
     var rangeLenRow = ((rhigh-rlow)/rinc).asInstanceOf[Int]+1    // row length of target range
     var rangeLenCol = ((chigh-clow)/cinc).asInstanceOf[Int]+1    // col length of target range

     if (mrN != rangeLenRow)  {  // improper sizes of the involved row subranges for assignment
         println("target vector row subrange and right-hand size vector lengths do not agree, i.e. mrN =  "+mrN+", rangeLenRow = "+rangeLenRow)
         this   // return the vector unaltered
     }

     if (mrM != rangeLenCol)  {  // improper sizes of the involved column subranges for assignment
         println("target vector column subrange and right-hand size vector lengths do not agree, i.e. mrM =  "+mrM+", rangeLenRow = "+rangeLenRow)
         this   // return the vector unaltered
     }

      if (rhigh >= N || chigh >= M)  {   // dynamically increase the size of the matrix when subassigning out of its range
          var newN = (Matrix.matResizeFactor*rhigh).asInstanceOf[Int]
          var newM = (Matrix.matResizeFactor*chigh).asInstanceOf[Int]
              // create the new larger matrix
          var v2 = new Array[Array[Double]](newN)
          var c=0
          while (c<newN) {
              v2(c) = new Array[Double](newM)
              c += 1
          }
            // copy the old one
            var i=0;  var j=0
            while (i<N)  {
                j=0;
                while (j<M)   {
                    v2(i)(j) = v(i)(j)
                    j += 1
                  }
                i += 1
            }
            // update the matrix representation to  the larger one
          v = v2
          N = newN
          M = newM
      }   // dynamically increase the size of the matrix


     // copy the values of the mr
        var rrow=1; var rcol=1; var lrowidx=rlow; var lcolidx = clow
        while (rrow < mr.N) {   // for all rows of the right-hand side matrix
            rcol=1
            lcolidx = clow   // starting column within the "subassigned" matrix
            while (rcol < mr.M)  {   // for all cols of the right-hand side matrix
                v(lrowidx)(lcolidx) = mr(rrow, rcol)
                lcolidx += cinc
                rcol += 1
            }
            lrowidx += rinc
            rrow += 1
          }

        }



    // Matrix -< Vec : column oriented addition
def -< (that: Vec): Matrix =  {
    var vl = that.length
    if (vl != this.N) 
      return this

      var nv = new Matrix(this.N, this.M)
   var i=1; var j=1
   while (j<=M) {   // for all columns
       i=1
    while (i<=N) {  // for all rows 
      nv.v(i)(j) = v(i)(j)-that(i-1)
      i +=1
    }
    j += 1
   }
   return nv
}

   // Matrix +< Vec : column oriented addition
def +< (that: Vec): Matrix =  {
    var vl = that.length
    if (vl != this.N) 
      return this

      var nv = new Matrix(this.N, this.M)
   var i=1; var j=1
   while (j<=M) {   // for all columns
       i=1
    while (i<=N) {  // for all rows 
      nv.v(i)(j) = v(i)(j)+that(i-1)
      i +=1
    }
    j += 1
   }
   return nv
}

   // Mat (NXM)*< Vec(MX1) : Mat(NX1) Matrix-Vector multiplication
def *< (that: Vec): Matrix =  {
    var vl = that.length
    if (vl != this.M ) 
      return this
   
    var sm=0.0
    var nv = new Matrix(this.N, 1)
   for (r<-1 to N) {      // all rows of the Mat
     sm=0.0
     for (c<-1 to M-1)
       sm += v(r)(c)*that(c-1)
     nv(r,1) = sm
    }
   
   nv
}

// dot product of a Matrix with a Matrix
def  dot(that: Matrix):  Matrix = {
  var nv = new Matrix(this.N-1, this.M-1 )
   var i=1; var j=1
   while (i<N) {
     j=1
    while (j<M)  {
      nv.v(i)(j) = v(i)(j) * that.v(i)(j)
      j += 1
    }
    i += 1
   }
 return nv
}

// dot product of a Matrix with an Array[Array[Double]]
def  dot(that: Array[Array[Double]]):  Matrix = {
  var nv = new Matrix(this.N-1, this.M-1 )
   var i=1; var j=1
   while (i<N) {
     j=1
    while (j<M)  {
      nv.v(i)(j) = v(i)(j) * that(i)(j)
      j += 1
    }
    i += 1
   }
 return nv
}


    // ROW - APPEND MATRIX TO MATRIX
// append rowwise the Matrix M
def  ::(rowsToAppend: Matrix): Matrix = {
    if (rowsToAppend.M != M )   // incompatible number of columns
      return this
    // create a new extended matrix to have also the added rows
    var  exrows = N+rowsToAppend.N-2   // new number of rows
    var  res = new Matrix(exrows, M-1)
    // copy "this" Matrix
    for (r<-1 to N-1)
       for (c<-1 to M-1)
         res.v(r)(c) = v(r)(c)
     println("copying first matrix")
    for (r<-0 to rowsToAppend.N-2)
       for (c<-1 to rowsToAppend.M-1)
         res.v(N+r)(c) = rowsToAppend.v(r+1)(c)
    res
}


    // COL - APPEND MATRIX TO MATRIX
// append columnwise the Matrix M
def  >::(colsToAppend: Matrix): Matrix = {
    if (colsToAppend.N != N )   // incompatible number of rows
      return this
    // create a new extended matrix to have also the added columns
    var  excols = M+colsToAppend.M-2   // new number of columns
    var  res = new Matrix(N-1, excols)
    // copy "this" Matrix
    for (r<-1 to N-1)
       for (c<-1 to M-1)
         res.v(r)(c) = v(r)(c)
    for (r<-1 to colsToAppend.N-1)
       for (c<-0 to colsToAppend.M-2)
         res.v(r)(M+c) = colsToAppend.v(r)(c+1)
    res
}


    // ROW - PREPEND MATRIX TO MATRIX
// prepend rowwise the Matrix M
def  :::(rowsToAppend: Matrix): Matrix = {
    if (rowsToAppend.M != M )   // incompatible number of columns
      return this
    // create a new extended matrix to have also the added rows
    var  exrows = N+rowsToAppend.N-2   // new number of rows
    var  res = new Matrix(exrows, M-1)
    // copy "this" Matrix 
    for (r<-1 to rowsToAppend.N-1)
       for (c<-1 to rowsToAppend.M-1)
         res.v(r)(c) = rowsToAppend.v(r)(c)
     println("copying rows ")
    for (r<-0 to N-2)
       for (c<-1 to M-1)
         res.v(rowsToAppend.N+r)(c) = v(r+1)(c)

    res
}

    // COL - PREPEND MATRIX TO MATRIX
// prepend columnwise the Matrix M
def  >:::(colsToAppend: Matrix): Matrix = {
    if (colsToAppend.N != N )   // incompatible number of rows
      return this
    // create a new extended matrix to have also the added columns
    var  excols = M+colsToAppend.M-2   // new number of columns
    var  res = new Matrix(N-1, excols)
    // copy "this" Matrix
    for (r<-1 to colsToAppend.N-1)
       for (c<-1 to colsToAppend.M-1)
         res.v(r)(c) = colsToAppend.v(r)(c)
    for (r<-1 to N-1)
       for (c<-0 to M-2)
         res.v(r)(colsToAppend.M+c) = v(r)(c+1) 

    res
}


// IN-PLACE Operations: Update directly the receiver, avoiding creating a new return object

    // Matrix + Matrix
def ++ (that: Matrix): Matrix =  {
  if (N != that.N || M != that.M)  // incompatible dimensions
      this
     else {
         var i=1; var j=1;
        while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) += that.v(i)(j)
      j +=1
    }
    i += 1
   }
   this
  }
}

    // Matrix + Double
def ++ (that: Double): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) += that
      j +=1
    }
    i += 1
   }
   this
}

 // Matrix + Int
def ++ (that: Int): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) += that
      j +=1
    }
    i += 1
   }
   this
}


    // Matrix - Matrix
def -- (that: Matrix): Matrix =  {
  if (N != that.N || M != that.M)  // incompatible dimensions
      this
     else {
         var i=1; var j=1;
        while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) -= that.v(i)(j)
      j +=1
    }
    i += 1
   }
   this
  }
}

    // Matrix - Double
def -- (that: Double): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) -= that
      j +=1
    }
    i += 1
   }
   this
}

 // Matrix - Int
def -- (that: Int): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) -= that
      j +=1
    }
    i += 1
   }
   this
}

     // Matrix - Long
def -- (that: Long): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) -= that
      j +=1
    }
    i += 1
   }
   this
}


    // Matrix * Double
def ** (that: Double): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) *= that
      j +=1
    }
    i += 1
   }
   this
}

    // Matrix * Int
def ** (that: Int): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) *= that
      j +=1
    }
    i += 1
   }
   this
}


    // Matrix * Long
def ** (that: Long): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) *= that
      j +=1
    }
    i += 1
   }
   this
}



    // Matrix / Double
def /| (that: Double): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) /= that
      j +=1
    }
    i += 1
   }
   this
}

    // Matrix / Int
def  /| (that: Int): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) /= that
      j +=1
    }
    i += 1
   }
   this
}


    // Matrix / Long
def /| (that: Long): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) /= that
      j +=1
    }
    i += 1
   }
   this
}




    // END OF IN-PLACE OPERATIONS
     // Matrix + Long
def ++ (that: Long): Matrix =  {
  var i=1; var j=1
  while (i<N) {
       j=1
    while (j<M) {
      v(i)(j) += that
      j +=1
    }
    i += 1
   }
   this
}



    // Matrix + Double
def + (that: Double): Matrix =  {
      var nv = new Matrix(this.N-1, this.M-1)
   var i=1; var j=1
   while (i<N) {
       j=1
    while (j<M) {
      nv.v(i)(j) = v(i)(j)+that
      j +=1
    }
    i += 1
   }
   return nv
}

  // Matrix - Double
def - (that: Double): Matrix =  {
      var nv = new Matrix(this.N-1, this.M-1)
   var i=1; var j=1
   while (i<N) {
       j=1
    while (j<M) {
      nv.v(i)(j) = v(i)(j)-that
      j +=1
    }
    i += 1
   }
   return nv
}

  
  // Matrix * Double
def * (that: Double): Matrix =  {
      var nv = new Matrix(this.N-1, this.M-1)
   var i=1; var j=1
   while (i<N) {
       j=1
    while (j<M) {
      nv.v(i)(j) = v(i)(j)*that
      j +=1
    }
    i += 1
   }
   return nv
}

    // Matrix + Double, for example at:   var m = ones(3,4); var mm = 10.0 +: m,
  // the implicit conversion of 10.0 to RichNumber is avoided since the expression is evaluated as m.+(10.0)
def +: (that: Double): Matrix =  {
      var nv = new Matrix(this.N-1, this.M-1)
   var i=1; var j=1
   while (i<N) {
       j=1
    while (j<M) {
      nv.v(i)(j) = v(i)(j)+that
      j +=1
    }
    i += 1
   }
   return nv
}

def -: (that: Double): Matrix =  {
      var nv = new Matrix(this.N-1, this.M-1)
   var i=1; var j=1
   while (i<N) {
       j=1
    while (j<M) {
      nv.v(i)(j) = v(i)(j)-that
      j +=1
    }
    i += 1
   }
   return nv
}



// Matrix * Double
def *: (that: Double): Matrix =  {
      var nv = new Matrix(this.N-1, this.M-1)
   var i=1; var j=1
   while (i<N) {
     j=1
    while (j<M) {
      nv.v(i)(j) = v(i)(j)*that
      j += 1
    }
    i += 1
   }
   return nv
}

    // unary Minus applied to a Matrix implies negation of all of its elements
def unary_- : Matrix =  {
      var nv = new Matrix(this.N-1, this.M-1)  // get a Matrix of the same dimension
   var i=1; var j=1
   while (i<N) {
     j=1
    while (j<M) {
      nv.v(i)(j) = -v(i)(j)  // negate element
      j += 1
    }
    i += 1
   }
   return nv
}

// transpose the Matrix
def trans: Matrix = {
   var nv = new Matrix(this.M-1, this.N-1)  // get a Matrix of dimension MXN
   var i=1; var j=1
   while (i<N) {
     j=1
    while (j<M) {
      nv.v(j)(i) = v(i)(j)  // negate element
      j += 1
    }
    i += 1
   }

   return nv
}


// transpose the Matrix
def T: Matrix = {
   var nv = new Matrix(this.M-1, this.N-1)  // get a Matrix of dimension MXN
   var i=1; var j=1;
   while (i<N) {
     j=1
    while (j<M) {
      nv.v(j)(i) = v(i)(j)  // negate element
      j += 1
    }
    i += 1
   }

   return nv
}

   

// Matrix + Matrix
def  + (that: Matrix): Matrix =  {
   var rN = this.N;   var rM = this.M;  // receiver's dimensions
   var sN = that.N;  var sM = that.M;  // parameter's dimensions

  if (rN == sN && rM == sM)  {     // same dimensions
       var nv = new Matrix(sN-1, sM-1)
       var i=1; var j=1;
   while (i<sN)  {
        j=1
       while (j<sM)  {
           nv.v(i)(j) = v(i)(j) + that.v(i)(j)
           j+=1
     }
     i += 1
   }
   return nv
  }
  else { // incompatible dimensions
   that
     }
}




// Matrix - Matrix
def  -  (that: Matrix): Matrix =  {
   var rN = this.N;   var rM = this.M;
   var sN = that.N;  var sM = that.M;
  
    if (rN == sN && rM == sM)  {     // same dimensions
       var nv = new Matrix(sN-1, sM-1)
       var i=1; var j=1;
   while (i<sN)  {
       j=1
    while (j<sM)  {
           nv.v(i)(j) = v(i)(j) - that.v(i)(j)
           j += 1
          }
         i += 1
    }
   return nv
  }  // same dimensions
  else { // incompatible dimensions
   return  that
     }
}

// Matrix * Matrix
 def * (that: Matrix): Matrix =  {
   var  rN = this.N;   var rM = this.M;
   var  sN = that.N;  var sM = that.M;
   
    var  v1Colj = new Array[Double](rM)
   var result = new Matrix(this.N-1, that.M-1)
   var j=1; var k=1;
   while (j < sM)  {
       k=1
      while  (k < rM) {
        v1Colj(k) =that.v(k)(j)
        k += 1
      }
    
      var i=1;
      while (i<rN) {
        var   Arowi = this.v(i)
        var   s = 0.0;
        k=1
        while (k< rM) {
          s += Arowi(k)*v1Colj(k)
          k += 1
        }
      result.v(i)(j) = s;
      i += 1
      }
 j += 1
   }
  return result
  }


// Matrix + Array[Array[Double]]
 def + (that: Array[Array[Double]]): Matrix =  {
   var  rN = this.N;   var rM = this.M;
   var  sN = that.length;  var sM = that(0).length

  if (rN == sN && rM == sM)  {     // same dimensions
       var nv = new Matrix(sN-1, sM-1)
       var i=1; var j=1;
   while (i<sN)  {
        j=1
       while (j<sM)  {
           nv.v(i)(j) = v(i)(j) + that(i)(j)
           j+=1
     }
     i += 1
   }
   return nv
  }
  else { // incompatible dimensions
   this
     }
}



// Matrix - Array[Array[Double]]
 def - (that: Array[Array[Double]]): Matrix =  {
   var  rN = this.N;   var rM = this.M;
   var  sN = that.length;  var sM = that(0).length

  if (rN == sN && rM == sM)  {     // same dimensions
       var nv = new Matrix(sN-1, sM-1)
       var i=1; var j=1
   while (i<sN)  {
        j=1
       while (j<sM)  {
           nv.v(i)(j) = v(i)(j) - that(i)(j)
           j+=1
     }
     i += 1
   }
   return nv
  }
  else { // incompatible dimensions
   this
     }
}


// Matrix * Array[Array[Double]]
 def * (that: Array[Array[Double]]): Matrix =  {
   var  rN = this.N;   var rM = this.M;
   var  sN = that.length;  var sM = that(0).length

  var  vr:Matrix =new Matrix(rN-1, sM-1)

  
   var   v1Colj = new Array[Double](rM)

   var j=1; var k=1;
   while (j < sM)  {
       k=1
      while  (k < rM) {
        v1Colj(k) =that(k)(j)
        k += 1
      }

      var i=1;
      while (i<rN) {
        var   Arowi = this.v(i)
        var   s = 0.0;
        k=1
        while (k< rM) {
          s += Arowi(k)*v1Colj(k)
          k += 1
        }
      vr.v(i)(j) = s;
      i += 1
      }
 j += 1
   }
  return vr
  }



//   the  M(fromR:incR:toR, fromC:incC:toC)  operation
def  subm(fromR:Int,  incR: Int,  toR: Int, fromC:Int,  incC: Int,  toC: Int)  =  
    apply(fromR, incR, toR, fromC, incC, toC)

//   the  M(fromR:incR:toR, :)  operation
def submr(fromR: Int, incR: Int, toR: Int) = 
    apply(':', fromR, incR, toR)

//   the  M(:, fromC:incC:toC)  operation
def submc(fromC: Int, incC: Int, toC: Int) = 
    apply(fromC, incC, toC, ':')


def ~() =  {
var  transposed = new Matrix(this.M-1, this.N-1)
   for (r<-1 to this.M-1)
    for (c<-1 to this.N-1)
      transposed(r,c) = this(c,r)
   transposed
} 


}


// Matrix's companion object
  object Matrix  {

   
 var mxElemsToDisplay = 6
 var matResizeFactor = 1.5
 

// construct a Matrix from a String 
// var m = M1("3.4 -6.7; -1.2 5.6")
def M1(s: String) =  {
    var nRows = 1
    var nCols = 0
    for (i<-0 to s.length-1)   // count how many rows are specified
      if  (s(i)==';')
        nRows += 1

  // seperate rows to an ArrayBuffer
   var buf = new scala.collection.mutable.ArrayBuffer[String]()
   var strtok = new java.util.StringTokenizer(s, ";")  // rows are separated by ';'
   while (strtok.hasMoreTokens) {
         val tok = strtok.nextToken
         buf +=tok
      }    

// count how many numbers each row has. Assuming that each line has the same number of elements
 val firstLine = buf(0)
 strtok = new java.util.StringTokenizer(firstLine, ", ")  // elements are separated by ',' or ' '
 while (strtok.hasMoreTokens) {
   val tok = strtok.nextToken
   nCols += 1  
}

// allocate double array
var elems = new Array[Array[Double]](nRows, nCols)
  // pass elements from ArrayBuffer to double array
for (k <- 0 to buf.length-1)  {  // read array
   var currentLine = buf(k)
   strtok = new java.util.StringTokenizer(currentLine, ", ")  // elements are separated by ',' or ' '
var c=0 
while (strtok.hasMoreTokens) {  // read row
   val tok = strtok.nextToken
   var cNumber = tok.toDouble
   elems(k)(c) = cNumber
   c += 1
     }   // read row
   }  // read array

 // return the converted Matrix
 new Matrix(elems)
 }  



//  Reshapes a matrix a to new dimension n X m
def  reshape(a: Matrix,  n: Int, m: Int): Matrix = {
   var  aCols  = a.M   // columns of matrix 
   var  aRows = a.N   // rows of matrix
   var N = n; var M = m;
  // m,n must be positive
 if (M<=0 || N<=0)  {  N = 1; M = aRows*aCols;  }  // default to reshaping in a large row
 if (N*M != aRows*aCols)  {
     return   a  // invalid new size: return the original matrix
    }
 var  nm = new Matrix(N, M)  // create the new matrix
// keep two set of indices, i.e. iorig, jorig: indices at the original matrix and inew, jnew: indices at the new matrix
var iorig=1; var jorig=1; var inew=1; var jnew=1
while (iorig < aRows)  {
    while (jorig < aCols)  {
        nm(inew, jnew) = a(iorig, jorig)
        jorig += 1
        jnew += 1
        if (jnew>=M) {  // next row of the reshaped new matrix
            jnew = 1
            inew += 1
           }
      } // jorig < aCols
    iorig += 1
    jorig = 1
  } // iorig < aRows
nm  // return the new Matrix
}


// resamples the Matrix every n rows and every m cols
def resample(matr: Matrix, n: Int,  m: Int):Matrix = {
    var matArr = matr.v
    var rows = matr.N
    var cols = matr.M
    var  rRows =  (rows/n).asInstanceOf[Int]-1
    var  rCols =  (cols/m).asInstanceOf[Int]-1
    var  newMatArr = new Array[Array[Double]](rRows, rCols) 
    var r=1; var c=1
    while (r<rRows) {
        c = 1
        while  (c<rCols) {
            newMatArr(r-1)(c-1) = matr.v(r*n)(m*c)
            c += 1
          }
      r += 1
}
   var   resampledMat = new Matrix(newMatArr);
    resampledMat
}

 // convert Vec to Matrix
def vecToMatrix( that: Vec): Matrix = {   
  var N = that.length
  var rmat = new Matrix(N, 1)
  for (k<-1 to N)
   rmat(k, 1) = that(k-1)

 rmat
}

 // COLUMN BASED ROUTINES
// columnwise sum
def sum(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var sm = 0.0
    var res = new Array[Double](M)
    for (ccol <-1  to M-1) {
     sm=0.0
     for (crow<- 1 to N-1)
       sm += matr.v(crow)(ccol)
     res(ccol) = sm 
     }
    res
}

// columnwise mean
def mean(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var sm = 0.0
    var res = new Array[Double](M)
    for (ccol <- 1  to M-1) {
     sm=0.0
     for (crow<- 1 to N-1)
       sm += matr.v(crow)(ccol)
     res(ccol) = sm/(N-1)
     }
    res
}

// columnwise product
def prod(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var pd = 1.0
    var res = new Array[Double](M)
    for (ccol <-1  to M-1) {
     pd=1.0
     for (crow<- 1 to N-1)
       pd *= matr.v(crow)(ccol)
     res(ccol) = pd 
     }
    res
}


// columnwise min
def min(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var res = new Array[Double](M)
    for (ccol <-1  to M-1) {
     var mn = matr.v(1)(ccol)  
     for (crow<- 2 to N-1)
        {
       var tmp = matr.v(crow)(ccol)
       if (tmp < mn)  mn = tmp
       }
     res(ccol) = mn
     }
    res
}


// columnwise max
def max(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var res = new Array[Double](M)
    for (ccol <-1  to M-1) {
     var mx = matr.v(1)(ccol)  
     for (crow<- 2 to N-1)
        {
       var tmp = matr.v(crow)(ccol)
       if (tmp > mx)  mx = tmp
       }
     res(ccol) = mx
     }
    res
}

// ROW BASED ROUTINES

// rowwise sum
def sumR(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var sm = 0.0
    var res = new Array[Double](N)
    for (crow<- 1 to N-1) {
     sm=0.0  // accumulates sum of all row elements
     for (ccol <-1  to M-1) 
        sm += matr.v(crow)(ccol)
     res(crow) = sm
      }
    res
}

// rowwise mean
def meanR(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var sm = 0.0
    var res = new Array[Double](N)
    for (crow<- 1 to N-1) {
     sm=0.0
    for (ccol <- 1  to M-1) 
       sm += matr.v(crow)(ccol)
     res(crow) = sm/(M-1)
     }
    res
}

// rowwise product
def prodR(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var pd = 1.0
    var res = new Array[Double](N)
    for (crow <-1  to N-1) {
     pd=1.0
     for (ccol<- 1 to M-1)
       pd *= matr.v(crow)(ccol)
     res(crow) = pd
     }
    res
}


// rowwise min
def minR(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var res = new Array[Double](N)
    for (crow <-1  to N-1) {
     var mn = matr.v(crow)(1)
     for (ccol<- 2 to M-1)
        {
       var tmp = matr.v(crow)(ccol)
       if (tmp < mn)  mn = tmp
       }
     res(crow) = mn
     }
    res
}


// rowwise max
def maxR(matr: Matrix): Array[Double] = {
    var N = matr.N;     var M = matr.M
    var res = new Array[Double](N)
    for (crow <-1  to N-1) {
     var mx = matr.v(crow)(1)
     for (ccol<- 2 to M-1)
        {
       var tmp = matr.v(crow)(ccol)
       if (tmp > mx)  mx = tmp
       }
     res(crow) = mx
     }
    res
}


  def  variance( v: Matrix): Array[Double] = {  // TODO: check
    var n = v.N; var m = v.M;
    var vvar = new Array[Double](m) // columns
    var degrees = n-1
    var c = 0.0; var s = 0.0
    var j=1
    while (j < m) { // all columns
        c = 0.0; s = 0;
        var k=1
        while (k<n) {  // all rows
            s += v(k,j)
            k += 1
           }
        s /= n
        k=1
        while (k<n)  {   // all rows
            c += ( v(k, j)-s)*(v(k, j)-s)
            k += 1
            }
        vvar(j) = c / degrees
        j += 1
     }   // all columns
    vvar
}


 def stddeviation( v: Matrix): Array[Double] = {  // TODO: check
   var varOfv = variance(v)
   var i=0
   while  (i<varOfv.length) {
     varOfv(i) = Math.sqrt(varOfv(i))
     i+=1
     }
 varOfv
    }

  
def  covariance( v1: Matrix, v2: Matrix):Matrix  = {   // TODO: check
    var n = v1.N
    var m1 = v1.M
    var m2 = v2.M
    var X = new Matrix(m1, m2)
    var degrees = n-1
    var s1=0.0; var s2 = 0.0; var c = 0.0;
    var i=1; var j=1
    while (i< m1)  {
        j=1
        while  (j<m2) {
            c = 0.0;   s1 = 0.0; s2 = 0.0
            var k=1
            while (k<n)  {
                s1 += v1(k, i)
                s2 += v2(k, i)
                k+=1
            }
          s1 /= m1;  s2 /= m2;
          k=1
          while  (k<n) {    c += (v1(k, i)-s1)*(v2(k, j)-s2); k+=1; }
         X(i, j) = c / degrees
          j+=1
        }
        i+=1
    }
    X    
    }



def  covariance( v: Matrix):Matrix  = {   // TODO: check
    var n = v.N
    var m = v.M
    var X = new Matrix(m, m)
    var degrees = n-1
    var s1=0.0; var s2 = 0.0; var c = 0.0;
    var i=1; var j=1
    while (i<m)  {
        j=1
        while  (j<m) {
            c = 0.0;   s1 = 0.0; s2 = 0.0
            var k=1
            while  (k<n)  {
                s1 += v(k, i)
                s2 += v(k, i)
                k += 1
            }
         s1 /= m;  s2 /= m;
         k=0
         while (k<n)  {
             c += (v(k, i)-s1)*(v(k, j)-s2)
             k += 1
         }
         X(i, j) = c / degrees
         j+=1
        }
        i+=1
     }
    X
    }

    

  def diag(n: Int): Matrix = {
         var om = new Matrix(n, n)
         var i=1;
         while  (i<= n) {
               om.v(i)(i) = 1.0
               i += 1
              }
          om
    }

def diag(a: Array[Double], oneIndexed: Boolean): Matrix = {
         var n = a.length-1
         var om = new Matrix(n, n)
         var i=1;
         while  (i<= n) {
               om.v(i)(i) = a(i)
               i += 1
              }
          om
    }

def diag(a: Array[Double]): Matrix = {
         var n = a.length
         var om = new Matrix(n, n)
         var i=1
         while  (i<= n) {
               om.v(i)(i) = a(i-1)
               i += 1
              }
          om
    }

     def ones(n: Int, m:Int): Matrix = {
         var om = new Matrix(n, m)
         var i=1; var j=1;
         while  (i <= n) {
            j=1
            while (j <= m) {
               om.v(i)(j) = 1.0
               j += 1
              }
            i += 1
         }
          om
    }

      def zeros(n: Int, m:Int): Matrix = {
         var om = new Matrix(n, m)
         var i=1; var j=1;
         while  (i<= n) {
            j=1
            while (j < m) {
               om.v(i)(j) = 0.0
               j += 1
              }
            i += 1
         }
          om
    }

 def fill(n: Int, m:Int, vl: Double): Matrix = {
         var om = new Matrix(n, m)
         var i=1; var j=1;
         while  (i<= n) {
            j=1
            while (j <= m) {
               om.v(i)(j) = vl
               j += 1
              }
            i += 1
         }
          om
    }

    def rand(n: Int, m:Int): Matrix = {
       println("IN rand() for Matrix")
       var om = new Matrix(n, m)
       var i=1; var j=1;
       while  (i<= n) {
            j=1
            while (j <= m) {
             om.v(i)(j) = java.lang.Math.random
             j += 1
           }
          i += 1
       }
       println("OUT OF rand() for Matrix")
       
          om
    }

  def randt(n:Int, m: Int): Matrix = {
    println("randj java")
    var matr = new Matrix(1,1)
    matr.N = n+1; matr.M = m+1
    var values = scalaSci.math.array.FMJ.randj(n+1, m+1)
    matr.v = values
    matr
  } 

    def randtest(n: Int, m:Int): Matrix = { //SOSSOS
       var values = new Array[Array[Double]](n+1, m+1)
       var i=1; var j=1;
       while  (i<= n) {
            j=1
            while (j <= m) {
             values(i)(j) = java.lang.Math.random
             j += 1
           }
          i += 1
       }
     var om = new Matrix(1, 1)
     om.setv(values, n, m)
     om
    }

def eye(n: Int): Matrix = {
   eye(n,n)
}

def  eye(n:Int, m:Int): Matrix = {
     var om = new Matrix(n, m)
     var minCoord = n
     if ( m<minCoord)
        minCoord = m
     var i=1 
     while  (i<= minCoord) {
             om.v(i)(i) = 1.0
              i += 1
         }
          om
     }


 def sin(v: Matrix): Matrix = {
       var N = v.N; var M = v.M;
        var om = new Matrix(N-1, M-1)
       var i=1; var j=1;
       while  (i< N) {
            j=1
            while (j < M) {
               om.v(i)(j) = java.lang.Math.sin(v.v(i)(j))
               j += 1
            }
            i += 1
       }
          om
    }


    def cos(v: Matrix): Matrix = {
        var N = v.N; var M = v.M;
        var om = new Matrix(N-1, M-1)
        var i=1; var j=1;
        while (i < N) {
            j=1
            while (j < M)  {
             om.v(i)(j) = java.lang.Math.cos(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

def tan(v: Matrix): Matrix = {
      var N = v.N; var M = v.M;
        var om = new Matrix(N-1, M-1)
        var i=1; var j=1;
        while (i < N) {
            j=1
            while (j < M)  {
             om.v(i)(j) = java.lang.Math.tan(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

    

def asin(v: Matrix): Matrix = {
      var N = v.N; var M = v.M;
        var om = new Matrix(N-1, M-1)
        var i=1; var j=1;
        while (i < N) {
            j=1
            while (j < M)  {
             om.v(i)(j) = java.lang.Math.asin(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }


def acos(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
           j=1
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.acos(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }




def atan(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
           j=1
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.atan(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

def sinh(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
           j=1
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.sinh(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



def cosh(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
          j=1
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.cosh(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



def tanh(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
          j=1
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.tanh(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
 }


def exp(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
           j=1
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.exp(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

def log(v: Matrix): Matrix = {
  var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
          j=1
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.log(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }


 def log2(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
var conv = java.lang.Math.log(2.0)
   var i=1; var j=1;
     while (i < N) {
           j=1
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.log(v.v(i)(j))/conv
             j += 1
            }
            i += 1
        }
          om
    }

def  log10(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
         j=1
         while (j < M)  {
             om.v(i)(j) = java.lang.Math.log10(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }


def abs(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
          j=1
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.abs(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



 def ceil(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
         j=1
         while (j < M)  {
             om.v(i)(j) = java.lang.Math.ceil(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }

  


def floor(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
          j=1
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.floor(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }




 def round(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
         j=1
         while (j < M)  {
             om.v(i)(j) = java.lang.Math.round(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }


   

def sqrt(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
          j=1
          while (j < M)  {
             om.v(i)(j) = java.lang.Math.sqrt(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }



 def toDegrees(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
           j=1
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.toDegrees(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }


 def toRadians(v: Matrix): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
         j=1
         while (j < M)  {
             om.v(i)(j) = java.lang.Math.toRadians(v.v(i)(j))
             j += 1
            }
            i += 1
        }
          om
    }




 def pow(v: Matrix, exponent:Double): Matrix = {
   var N = v.N; var M = v.M;
   var om = new Matrix(N-1, M-1)
   var i=1; var j=1;
     while (i < N) {
           j=1
           while (j < M)  {
             om.v(i)(j) = java.lang.Math.pow(v.v(i)(j), exponent)
             j += 1
            }
            i += 1
        }
          om
    }



def  transpose(Mt: Matrix): Matrix =  {
   var  transposed = new Matrix(Mt.M-1, Mt.N-1)
   for (r<-1 to Mt.M-1)
    for (c<-1 to Mt.N-1)
      transposed(r,c) = Mt(c,r)
   transposed
}

  def  T(Mt: Matrix): Matrix =  {
    transpose(Mt)
	}


def  LU(a: Matrix): (Matrix, Matrix, Matrix) = {
   var aux  = new Array[Double](4)
   var n = a.length
   var p = new Array[Int](n+1) 
   var P = new Matrix(a.N-1, a.M-1)
   var aa = new Array[Array[Double]](a.N, a.M)
   for (r<-1 until a.N)
    for (c<-1 until a.M)
      aa(r)(c) = a(r,c)

   aux(2) = 1.0e-5
  numal.Linear_algebra.dec(aa, n,  aux, p)

 for (k<-1 to n)
   println("p("+k+") == "+p(k))
 println()
 // construct permutation matrix
  for (k<-1 to n)  {
      P(k, k) = 1
      if (k!=p(k)) {
         P(k, k) = 0
         P(p(k), p(k)) = 0
         P(k, p(k)) = 1
         P(p(k), k) = 1
          }
        }
    
// extract lower and upper matrices from aa
  var lm = new  Array[Array[Double]](a.N, a.M)
  var um = new  Array[Array[Double]](a.N, a.M)
  for (r<-1 until a.N)
    for (c<-1 to r)
      lm(r)(c) = aa(r)(c)
  for (r<-1 until a.N) {
    um(r)(r)=1.0
    for (c<-r+1 until a.M )
      um(r)(c) = aa(r)(c)
   um(r)(r) = 1
  }


 val oneIndexed = true
 (new Matrix(lm, oneIndexed), new Matrix(um, oneIndexed), P)
}

 
  def det(a: Matrix): Double = {
    numal.Linear_algebra.det(a.v)
  }
// solves a well-conditioned linear system of equations Ax=b whose order is small
// relative to the number of binary digits in the number representation
def decsol( a: Matrix, aux: Array[Double], b: Array[Double]): Array[Double] =  {
  numal.Linear_algebra.decsol(a.v, a.N-1, aux, b)
  b
 }

    // copy in place the inverse matrix
 def decinv( a: Matrix): Matrix = {
     var aux = new Array[Double](4)
     aux(2) = 1.0e-5
     var ainv = new  Matrix(a)
     numal.Linear_algebra.decinv(ainv.v, ainv.N-1, aux)
     ainv
   }


 def decinv( a: Matrix, aux: Array[Double]): Matrix = {
     numal.Linear_algebra.decinv(a.v, a.N-1, aux)
     a
 }

def gsssol( a: Matrix, b: Array[Double]): Matrix = {
    var aux = new Array[Double](8)
    aux(2) = 1.0e-5
    aux(4) = 8
    numal.Linear_algebra.gsssol(a.v, a.N-1, aux, b)
    new Matrix(b)
}

def gsssol( a: Matrix, aux: Array[Double], b: Array[Double]): Array[Double] = {
    numal.Linear_algebra.gsssol(a.v, a.N-1, aux, b)
    b
}

def eigcom(ar:  Matrix, ai: Matrix, em: Array[Double])= {
    // : (valr: Array[Double], vali: Array[Double], vr: Array[Array[Double]], vi: Array[Array[Double]]) 
    var N = ar.N
    var dvalr = new Array[Double](N); var dvali = new Array[Double](N);
    var dvr = new Array[Array[Double]](N,N); var dvi = new Array[Array[Double]](N,N);

     numal.Linear_algebra.eigcom(ar.v, ai.v, N-1, em, dvalr, dvali, dvr, dvi )
     // return the results as Matrix'es for more convenient handling
     var mdvalr = new Matrix(dvalr)
     var mdvali = new Matrix(dvali)
     var mdvr = new Matrix(dvr)
     var mdvi = new Matrix(dvi)
        (mdvr, mdvi, mdvalr, mdvali)

}

  def eigcom(ar:  Matrix, ai: Matrix)= {
    // : (valr: Array[Double], vali: Array[Double], vr: Array[Array[Double]], vi: Array[Array[Double]])
    var N = ar.N
    var  em =  new Array[Double](8)
    em(0) = 5.0e-6;  em(2)=1.0e-5; em(4)=10.0; em(6)=10.0;
    var dvalr = new Array[Double](N); var dvali = new Array[Double](N+1);
    var dvr = new Array[Array[Double]](N,N); var dvi = new Array[Array[Double]](N,N);

     numal.Linear_algebra.eigcom(ar.v, ai.v, N-1, em, dvalr, dvali, dvr, dvi )
     var mdvalr = new Matrix(dvalr)
     var mdvali = new Matrix(dvali)
     var mdvr = new Matrix(dvr)
     var mdvi = new Matrix(dvi)
        (mdvr, mdvi, mdvalr, mdvali)
}


  def eigcom(ar:  Matrix)= {
    // : (valr: Array[Double], vali: Array[Double], vr: Array[Array[Double]], vi: Array[Array[Double]])
    var N = ar.N
    println(" N = "+N)
    var ai = zeros(N, N)
    var  em =  new Array[Double](8)
    em(0) = 5.0e-6;  em(2)=1.0e-5; em(4)=10.0; em(6)=10.0;
    var dvalr = new Array[Double](N); var dvali = new Array[Double](N);
    var dvr = new Array[Array[Double]](N,N); var dvi = new Array[Array[Double]](N,N);

     numal.Linear_algebra.eigcom(ar.v, ai.v, N-1, em, dvalr, dvali, dvr, dvi )
     var mdvalr = new Matrix(dvalr)
     var mdvali = new Matrix(dvali)
     var mdvr = new Matrix(dvr)
     var mdvi = new Matrix(dvi)
        (mdvr, mdvi, mdvalr, mdvali)
}

    def eig(ar:  Matrix, ai: Matrix, em: Array[Double])= {
    var N = ar.N
    var dvalr = new Array[Double](N); var dvali = new Array[Double](N);
    var dvr = new Array[Array[Double]](N,N); var dvi = new Array[Array[Double]](N,N);

     numal.Linear_algebra.eigcom(ar.v, ai.v, N-1, em, dvalr, dvali, dvr, dvi )
     // return the results as Matrix'es for more convenient handling
     var mdvalr = new Matrix(dvalr)
     var mdvali = new Matrix(dvali)
     var mdvr = new Matrix(dvr)
     var mdvi = new Matrix(dvi)
        (mdvalr, mdvali, mdvr, mdvi)

}

  def eig(ar:  Matrix, ai: Matrix)= {
    var N = ar.N
    var  em =  new Array[Double](8)
    em(0) = 5.0e-6;  em(2)=1.0e-5; em(4)=10.0; em(6)=10.0;
    var dvalr = new Array[Double](N); var dvali = new Array[Double](N+1);
    var dvr = new Array[Array[Double]](N,N); var dvi = new Array[Array[Double]](N,N);

     numal.Linear_algebra.eigcom(ar.v, ai.v, N-1, em, dvalr, dvali, dvr, dvi )
     var mdvalr = new Matrix(dvalr)
     var mdvali = new Matrix(dvali)
     var mdvr = new Matrix(dvr)
     var mdvi = new Matrix(dvi)
        (mdvalr, mdvali, mdvr, mdvi)
}


  def eig(ar:  Matrix)= {
    var N = ar.N-1
    println(" N = "+N)
    var arcp = new Array[Array[Double]](N+1, N+1)
      // make a copy of the matrix to pass to comeig1
    for (r<-1 to N)
      for (c<-1 to N) 
         arcp(r)(c) = ar.v(r)(c)

    var  em =  new Array[Double](10)
    em(0) = 5.0e-6;   // the machine precision
    em(2)=1.0e-5;   // the relative tolerance used for the QR iteration (em(2) > em(0))
    em(4)=10.0;   // the maximum allowed number of iterations
    em(6)=1.0e-5; // the tolerance used for the eigenvectors
    println("EM(6) ="+em(6))
    em(8)=5; // the maximum allowed number of inverse iterations for the calculation of each calculation
      // arrays for the real and imaginary parts of the calculated eigenvalues of the given matrix
    var re = new Array[Double](N+1); 
    var im = new Array[Double](N+1);
      // the calculated eigenvectors are delivered in the columns of vec; an eigenvector corresponding to a real eigenvalue given in array re
      // is delivered in the corresponding column of array vec;
      // the real and imaginary part of an eigenvector corresponding to the first member of a nonreal complex conjugate pair
     // of eigenvalues given in the arrays re, im are delivered in the two consecutive columns of array vec corresponding to this pair
     // (the eigenvectors corresponding to the second members of nonreal complex complex conjugate pairs are not delivered,
     // since they are simply the complex conjugate of those corresponding to the first member of such pairs)
    var vec = new Array[Array[Double]](N+1,N+1); 
    numal.Linear_algebra.comeig1(arcp, N, em, re, im, vec )

    var vecEig= new Array[Array[Double]](N, N)  // the eigenvectors
        // the calculated eigenvectors are delivered in the columns of vec
        // one-indexed is also the [][] vec result !!
       var r = 0; var c  = 0
     while (r< N) {
       c = 0
       while  (c< N) {
          vecEig(r)(c) = vec(r+1)(c+1)  // copy eigenvectors
          c += 1
       }
       r += 1
     }

    var vecEiv= new Array[Array[Double]](N, 2)  // the eigenvalues
    r = 0
    while  (r<N)  {
        vecEiv(r)(0) = re(r+1)  // real eigenvalues
        vecEiv(r)(1) = im(r+1)  // imaginary eigenvalues
        r += 1
      }
 (new Matrix(vecEig), new Matrix(vecEiv)  )

}

// computes the SVD of a given matrix A
def svdRank(A: Matrix) = {
  var M = A.Nrows
  var N = A.Ncols
  var singValues = new Array[Double](N+1)    // the singular values
  var U = new Array[Array[Double]](M+1, N+1)
  var V = new Array[Array[Double]](N+1, N+1)

// copy data matrix to U
for (r<-1 to M)
  for (c<-1 to N)
    U(r)(c) = A.v(r)(c)

var em = new Array[Double](8)
em(0) = 1.0e-6  // machine's precision
em(2) = 1.0e-5  // the relative precision in the singular values
em(4) = 25.0  // the maximal number of iterations to be performed
em(5) = 1.0e-5   // the minimal nonneglectable singular value

// call the NUMAL routine to perform the SVD
 numal.Linear_algebra.qrisngvaldec(U, M, N, singValues, V, em)

// Return
// em(1): the infinity norm of the matrix
// em(3): the maximal neglected superdiagonal element
// em(5): the number of iterations performed
// em(7): the numerical rank of the matrix; i.e. the number of singular values greater than or equal to em(6)

(new Matrix(U, true), diag(singValues, true), new Matrix(V, true), em(7), em(1), em(3), em(5))
}

// computes the SVD of a given matrix A
def svd(A: Matrix) = {
  
var S  = svdRank(A)
  (S._1, S._2, S._3)
}



 

// construct a Matrix from a double [][] array, zero indexed
def MatrixFromZeroIndexedDArr(vals: Array[Array[Double]]): Matrix = {
     new Matrix(vals)
 }


def testMatrix(N: Int, M: Int) = {
  var a = new Matrix(N, M)
  for (rows <- 1 to N)
    for (cols <-1 to M)
      a(rows, cols) = (rows-1)*10+cols-1

  a
}


def norm2(M: Matrix): Double = {
  LinearAlgebra.norm2( M.v)
}

def  norm1(M: Matrix): Double = {
  LinearAlgebra.norm1( M.v)
}

def  normF(M: Matrix): Double = {
  LinearAlgebra.normF( M.v)
}

def  normInf(M: Matrix): Double = {
   LinearAlgebra.normInf( M.v)
}
  

/*
var m = rand(5,5)
var (dvalr, dvali, dvr, dvi)= eigcom(m)
*/
}


